{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regression\n",
    "\n",
    "Different from classification, the output variable takes continuous values in regression analysis. Smile's regression algorithms are in the package smile.regression and all algorithms implement the interface Regression that has a single method predict to apply the model to an instance. For all algorithms, the model can be trained through the constructor. Meanwhile, each algorithm has a Trainer companion class that can hold model hyperparameters and be applied to multiple training datasets.\n",
    "\n",
    "The high-level operators are defined in Scala trait smile.regression.Operators and also in the package object of smile.regression. In what follows, we will discuss each algorithm, their high-level Scala API, and examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import $ivy.`com.github.haifengl::smile-scala:4.0.0`\n",
    "import $ivy.`org.slf4j:slf4j-simple:2.0.16`  \n",
    "\n",
    "import scala.language.postfixOps\n",
    "import smile._\n",
    "import smile.math._\n",
    "import smile.math.distance._\n",
    "import smile.math.kernel._\n",
    "import smile.math.matrix._\n",
    "import smile.math.matrix.Matrix._\n",
    "import smile.math.rbf._\n",
    "import smile.stat.distribution._\n",
    "import smile.data._\n",
    "import smile.data.formula._\n",
    "import smile.data.measure._\n",
    "import smile.data.`type`._\n",
    "import smile.base.cart.SplitRule\n",
    "import smile.base.mlp._\n",
    "import smile.base.rbf.RBF\n",
    "import smile.regression._\n",
    "import smile.feature._\n",
    "import smile.validation._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ordinary Least Squares\n",
    "\n",
    "In linear regression, the model specification is that the dependent variable is a linear combination of the parameters. The residual is the difference between the value of the dependent variable predicted by the model, and the true value of the dependent variable. Ordinary least squares obtains parameter estimates that minimize the sum of squared residuals, SSE (also denoted RSS).\n",
    "\n",
    "The ordinary least squares (OLS) estimator is consistent when the independent variables are exogenous and there is no multicollinearity, and optimal in the class of linear unbiased estimators when the errors are homoscedastic and serially uncorrelated. Under these conditions, the method of OLS provides minimum-variance mean-unbiased estimation when the errors have finite variances.\n",
    "```\n",
    "def ols(formula: Formula, data: DataFrame, method: String = \"qr\", stderr: Boolean = true, recursive: Boolean = true): LinearModel\n",
    "``` \n",
    "There are several different frameworks in which the linear regression model can be cast in order to make the OLS technique applicable. Each of these settings produces the same formulas and same results, the only difference is the interpretation and the assumptions which have to be imposed in order for the method to give meaningful results. The choice of the applicable framework depends mostly on the nature of data at hand, and on the inference task which has to be performed.\n",
    "\n",
    "Least squares corresponds to the maximum likelihood criterion if the experimental errors have a normal distribution and can also be derived as a method of moments estimator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "val planes = read.arff(\"../data/weka/regression/2dplanes.arff\")\n",
    "val model = lm(\"y\" ~ \".\", planes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then can apply the regression model on a new data point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.predict(planes(0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once a regression model has been constructed, it may be important to confirm the goodness of fit of the model and the statistical significance of the estimated parameters. Commonly used checks of goodness of fit include the R-squared, analysis of the pattern of residuals and hypothesis testing. Statistical significance can be checked by an F-test of the overall fit, followed by t-tests of individual parameters.\n",
    "\n",
    "Simply print out the model, we can inspect R-squared, the statistics and p-values of the t-test of parameter significance, and the F-test of goodness-of-fit. Interpretations of these diagnostic tests rest heavily on the model assumptions. Although examination of the residuals can be used to invalidate a model, the results of a t-test or F-test are sometimes more difficult to interpret if the model's assumptions are violated. For example, if the error term does not have a normal distribution, in small samples the estimated parameters will not follow normal distributions and complicate inference. With relatively large samples, however, a central limit theorem can be invoked such that hypothesis testing may proceed using asymptotic approximations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ridge Regression\n",
    "\n",
    "Coefficient estimates for multiple linear regression models rely on the independence of the model terms. When terms are correlated and the columns of the design matrix `X` have an approximate linear dependence, the matrix `X'X` becomes close to singular. As a result, the least-squares estimate becomes highly sensitive to random errors in the observed response Y, producing a large variance.\n",
    "\n",
    "Ridge regression is one method to address these issues. In ridge regression, the matrix `X'X` is perturbed so as to make its determinant appreciably different from 0.\n",
    "\n",
    "Ridge regression is a kind of Tikhonov regularization, which is the most commonly used method of regularization of ill-posed problems. Ridge regression shrinks the regression coefficients by imposing a penalty on their size. By allowing a small amount of bias in the estimates, more reasonable coefficients may often be obtained. Often, small amounts of bias lead to dramatic reductions in the variance of the estimated model coefficients.\n",
    "\n",
    "Another interpretation of ridge regression is available through Bayesian estimation. In this setting the belief that weight should be small is coded into a prior distribution.\n",
    "```\n",
    "def ridge(formula: Formula, data: DataFrame, lambda: Double): LinearModel\n",
    "``` \n",
    "where the parameter x is a matrix containing the explanatory variables, y is the response values, and lambda is the shrinkage/regularization parameter. Large lambda means more shrinkage. Choosing an appropriate value of lambda is important, and also difficult.\n",
    "\n",
    "Longley's macroeconomic regression data has 7 economical variables, observed yearly from 1947 to 1962. It is a well-known example for a highly collinear regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "val longley = read.arff(\"../data/weka/regression/longley.arff\")\n",
    "val model = ridge(\"employed\" ~ \".\", longley, 0.0057)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lasso Regression\n",
    "\n",
    "Lasso (least absolute shrinkage and selection operator) regression is a shrinkage and selection method for linear regression. It minimizes the usual sum of squared errors, with a bound on the sum of the absolute values of the coefficients (i.e. L1-regularized). It has connections to soft-thresholding of wavelet coefficients, forward stage-wise regression, and boosting methods.\n",
    "\n",
    "The Lasso typically yields a sparse solution, of which the parameter vector β has relatively few nonzero coefficients. In contrast, the solution of L2-regularized least squares (i.e. ridge regression) typically has all coefficients nonzero. Because it effectively reduces the number of variables, the Lasso is useful in some contexts.\n",
    "\n",
    "There is no analytic formula or expression for the optimal solution to the L1-regularized least squares problems. Therefore, its solution must be computed numerically. The objective function in the L1-regularized least squares is convex but not differentiable, so solving it is more of a computational challenge than solving the L2-regularized least squares. The Lasso may be solved using quadratic programming or more general convex optimization methods, as well as by specific algorithms such as the least angle regression algorithm.\n",
    "```\n",
    "def lasso(formula: Formula, data: DataFrame, lambda: Double, tol: Double = 1E-3, maxIter: Int = 5000): LinearModel\n",
    "``` \n",
    "where the parameter x is a matrix containing the explanatory variables, y is the response values, lambda is the shrinkage/regularization parameter, tol is the tolerance for stopping iterations (relative target duality gap), and maxIter is the maximum number of iterations.\n",
    "\n",
    "For over-determined systems (more instances than variables, commonly in machine learning), we normalize variables with mean 0 and standard deviation 1. For under-determined systems (less instances than variables, e.g. compressed sensing), we assume white noise (i.e. no intercept in the linear model) and do not perform normalization. Note that the solution is not unique in this case.\n",
    "\n",
    "In what follows, we will apply Lasso to the diabetes data used in the \"Least Angle Regression\" paper. The basic data has 10 columns standardized to have unit L2 norm in each column and zero mean. The data used in the below example has 64 columns that consists of basic data plus certain interactions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "val diabetes = read.csv(\"../data/regression/diabetes.csv\")\n",
    "val model = lasso(\"y\" ~ \".\", diabetes, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Radial Basis Function Networks\n",
    "\n",
    "A radial basis function network is an artificial neural network that uses radial basis functions as activation functions. It is a linear combination of radial basis functions. They are used in function approximation, time series prediction, and control.\n",
    "\n",
    "In its basic form, radial basis function network is in the form\n",
    "```\n",
    "    y(x) = Σ wi φ(||x-ci||)\n",
    "``` \n",
    "where the approximating function `y(x)` is represented as a sum of `N` radial basis functions `φ`, each associated with a different center `ci`, and weighted by an appropriate coefficient wi. For distance, one usually chooses Euclidean distance. The weights wi can be estimated using the matrix methods of linear least squares, because the approximating function is linear in the weights.\n",
    "\n",
    "RBF networks are typically trained by a two-step algorithm. In the first step, the center vectors ci of the RBF functions in the hidden layer are chosen. This step can be performed in several ways; centers can be randomly sampled from some set of examples, or they can be determined using k-means clustering. Note that this step is unsupervised.\n",
    "\n",
    "The second step simply fits a linear model with coefficients wi to the hidden layer's outputs with respect to some objective function. A common objective function, at least for regression/function estimation, is the least squares function.\n",
    "\n",
    "A optional third backpropagation step can be performed to fine-tune all of the RBF network's parameters.\n",
    "```\n",
    "def rbfnet[T](x: Array[T], y: Array[Double], neurons: Array[RBF[T]], normalized: Boolean): RBFNetwork[T]\n",
    "\n",
    "def rbfnet(x: Array[Array[Double]], y: Array[Double], k: Int, normalized: Boolean = false): RBFNetwork[Array[Double]]\n",
    "``` \n",
    "The popular choices for `φ` comprise the Gaussian function and the so called thin plate splines. The advantage of the thin plate splines is that their conditioning is invariant under scalings. Gaussian, multi-quadric and inverse multi-quadric are infinitely smooth and and involve a scale or shape parameter, `r0 > 0`. Decreasing `r0` tends to flatten the basis function. For a given function, the quality of approximation may strongly depend on this parameter. In particular, increasing `r0` has the effect of better conditioning (the separation distance of the scaled points increases)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "val y = diabetes(\"y\").toDoubleArray()\n",
    "val x = diabetes.select(1 until 11).toArray() // use only the primary attributes\n",
    "cv.regression(10, x, y) { case (x, y) => rbfnet(x, y, 10) }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Support Vector Regression\n",
    "\n",
    "Support vector machine can be used as a regression method, maintaining all the main features of the algorithm. In the case of regression, a margin of tolerance ε is set in approximation. The goal of SVR is to find a function that has at most ε deviation from the response variable for all the training data, and at the same time is as flat as possible. In other words, we do not care about errors as long as they are less than ε, but will not accept any deviation larger than this.\n",
    "\n",
    "Like SVM for classification, nonlinear SVR employs kernel trick for implict mapping. And the model produced by SVR depends only on a subset of the training data, because the cost function ignores any training data close to the model prediction (within the ε threshold).\n",
    "```\n",
    "def svr[T](x: Array[T],\n",
    "           y: Array[Double],\n",
    "           kernel: MercerKernel[T],\n",
    "           eps: Double,\n",
    "           C: Double,\n",
    "           tol: Double = 1E-3): KernelMachine[T]\n",
    "``` \n",
    "where the parameter `x` is the training data, y is the response variable, kernel is the kernel function, eps is the loss function error threshold, `C` is the soft margin penalty parameter, weight is the optional positive instance weight so that the soft margin penalty parameter for instance `i` will be `weight(i) * C`, and tol the tolerance of convergence test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv.regression(10, x, y) { case (x, y) => svm(x, y, new GaussianKernel(0.06), 20, 10) }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regression Tree\n",
    "\n",
    "Similar to decision trees, regression trees can be learned by splitting the training set into subsets based on an attribute value test. This process is repeated on each derived subset in a recursive manner called recursive partitioning. The recursion is completed when the subset at a node all has the same value of the target variable, or when splitting no longer adds value to the predictions.\n",
    "\n",
    "The algorithms that are used for constructing decision trees usually work top-down by choosing a variable at each step that is the next best variable to use in splitting the set of items. \"Best\" is defined by how well the variable splits the set into homogeneous subsets that have the same value of the target variable. Different algorithms use different formulae for measuring \"best\". Used by the CART algorithm, Gini impurity is a measure of how often a randomly chosen element from the set would be incorrectly labeled if it were randomly labeled according to the distribution of labels in the subset. Gini impurity can be computed by summing the probability of each item being chosen times the probability of a mistake in categorizing that item. It reaches its minimum (zero) when all cases in the node fall into a single target category. Information gain is another popular measure, used by the ID3, C4.5 and C5.0 algorithms. Information gain is based on the concept of entropy used in information theory. For categorical variables with different number of levels, however, information gain are biased in favor of those attributes with more levels. Instead, one may employ the information gain ratio, which solves the drawback of information gain.\n",
    "```\n",
    "def cart(formula: Formula,\n",
    "         data: DataFrame,\n",
    "         maxDepth: Int = 20,\n",
    "         maxNodes: Int = 0,\n",
    "         nodeSize: Int = 5): RegressionTree\n",
    "``` \n",
    "where the parameter x is the training data, y is the response variable, maxNodes is the maximum number of leaf nodes in the tree as a regularization, and The optional attributes is the attribute properties. If not provided, all attributes are treated as numeric values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv.regression(10, \"y\" ~ \".\", diabetes) { (formula, data) => cart(formula, data, 200) }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random Forest\n",
    "\n",
    "Random forest is an ensemble classifier that consists of many decision trees and outputs the majority vote of individual trees. The method combines bagging idea and the random selection of features.\n",
    "\n",
    "Each tree is constructed using the following algorithm:\n",
    "\n",
    "- If the number of cases in the training set is N, randomly sample N cases with replacement from the original data. This sample will be the training set for growing the tree.\n",
    "- If there are M input variables, a number m << M is specified such that at each node, m variables are selected at random out of the M and the best split on these m is used to split the node. The value of m is held constant during the forest growing.\n",
    "- Each tree is grown to the largest extent possible. There is no pruning.\n",
    "\n",
    "```\n",
    "def randomForest(formula: Formula,\n",
    "                 data: DataFrame,\n",
    "                 ntrees: Int = 500,\n",
    "                 mtry: Int = 0,\n",
    "                 maxDepth: Int = 20,\n",
    "                 maxNodes: Int = 0,\n",
    "                 nodeSize: Int = 5,\n",
    "                 subsample: Double = 1.0): RandomForest\n",
    "``` \n",
    "where ntrees is the number of trees, nodeSize is the number of instances in a node below which the tree will not split, setting `nodeSize = 5` generally gives good results, and maxNodes is the maximum number of leaf nodes. The parameter mtry is the number of input variables to be used to determine the decision at a node of the tree. The default value `p/3` seems to give generally good performance, where p is the number of variables. Besides, subsample is the sampling rate for training tree. The default value 1.0 means sampling with replacement. Otherwise, the value < 1.0 means sampling without replacement. The optional attributes is the attribute properties. If not provided, all attributes are treated as numeric values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv.regression(10, \"y\" ~ \".\", diabetes) { (formula, data) => randomForest(formula, data, maxNodes = 100) }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient Boosting\n",
    "\n",
    "Gradient boosting is typically used with decision trees (especially CART regression trees) of a fixed size as base learners. For this special case Friedman proposes a modification to gradient boosting method which improves the quality of fit of each base learner.\n",
    "\n",
    "Generic gradient boosting at the t-th step would fit a regression tree to pseudo-residuals. Let J be the number of its leaves. The tree partitions the input space into J disjoint regions and predicts a constant value in each region. The parameter J controls the maximum allowed level of interaction between variables in the model. With `J = 2` (decision stumps), no interaction between variables is allowed. With `J = 3` the model may include effects of the interaction between up to two variables, and so on. Hastie et al. comment that typically `4 ≤ J ≤ 8` work well for boosting and results are fairly insensitive to the choice of in this range, `J = 2` is insufficient for many applications, and J > 10 is unlikely to be required.\n",
    "\n",
    "Fitting the training set too closely can lead to degradation of the model's generalization ability. Several so-called regularization techniques reduce this over-fitting effect by constraining the fitting procedure. One natural regularization parameter is the number of gradient boosting iterations T (i.e. the number of trees in the model when the base learner is a decision tree). Increasing T reduces the error on training set, but setting it too high may lead to over-fitting. An optimal value of T is often selected by monitoring prediction error on a separate validation data set.\n",
    "\n",
    "Another regularization approach is the shrinkage which times a parameter η (called the \"learning rate\") to update term. Empirically it has been found that using small learning rates (such as η < 0.1) yields dramatic improvements in model's generalization ability over gradient boosting without shrinking (`η = 1`). However, it comes at the price of increasing computational time both during training and prediction: lower learning rate requires more iterations.\n",
    "```\n",
    "def gbm(formula: Formula,\n",
    "        data: DataFrame,\n",
    "        loss: Loss = Loss.lad(),\n",
    "        ntrees: Int = 500,\n",
    "        maxDepth: Int = 20,\n",
    "        maxNodes: Int = 6,\n",
    "        nodeSize: Int = 5,\n",
    "        shrinkage: Double = 0.05,\n",
    "        subsample: Double = 0.7): GradientTreeBoost\n",
    "``` \n",
    "where the parameter x is the training data, y is the response variable, ntrees is the number of trees, and maxNodes is the maximum number of leaf nodes. The parameter loss is the loss function for regression. By default, least absolute deviation is employed for robust regression. The parameter shrinkage is the shrinkage parameter in (0, 1] controls the learning rate of procedure. Besides, subsample is the sampling rate for training tree. The default value 1.0 means sampling with replacement. Otherwise, the value < 1.0 means sampling without replacement. The optional attributes is the attribute properties. If not provided, all attributes are treated as numeric values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv.regression(10, \"y\" ~ \".\", diabetes) { (formula, data) => gbm(formula, data) }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gaussian Process\n",
    "\n",
    "A Gaussian process is a stochastic process whose realizations consist of random values associated with every point in a range of times (or of space) such that each such random variable has a normal distribution. Moreover, every finite collection of those random variables has a multivariate normal distribution.\n",
    "\n",
    "A Gaussian process can be used as a prior probability distribution over functions in Bayesian inference. Given any set of N points in the desired domain of your functions, take a multivariate Gaussian whose covariance matrix parameter is the Gram matrix of N points with some desired kernel, and sample from that Gaussian. Inference of continuous values with a Gaussian process prior is known as Gaussian process regression.\n",
    "\n",
    "The fitting is performed in the reproducing kernel Hilbert space with the \"kernel trick\". The loss function is squared-error. This also arises as the kriging estimate of a Gaussian random field in spatial statistics.\n",
    "\n",
    "A significant problem with Gaussian process prediction is that it typically scales as `O(n3)`. For large problems (e.g. `n > 10,000`) both storing the Gram matrix and solving the associated linear systems are prohibitive on modern workstations. An extensive range of proposals have been suggested to deal with this problem. A popular approach is the reduced-rank Approximations of the Gram Matrix, known as Nystrom approximation. Greedy approximation is another popular approach that uses an active set of training points of size m selected from the training set of size n > m. We assume that it is impossible to search for the optimal subset of size m due to combinatorics. The points in the active set could be selected randomly, but in general we might expect better performance if the points are selected greedily w.r.t. some criterion. Recently, researchers had proposed relaxing the constraint that the inducing variables must be a subset of training/test cases, turning the discrete selection problem into one of continuous optimization.\n",
    "```\n",
    "object gpr {\n",
    "  // regular Gaussian process model\n",
    "  def apply[T](x: Array[T], y: Array[Double], kernel: MercerKernel[T], lambda: Double): KernelMachine[T]\n",
    "  // approximate Gaussian process model\n",
    "  def approx[T](x: Array[T], y: Array[Double], t: Array[T], kernel: MercerKernel[T], lambda: Double): KernelMachine[T]\n",
    "  // Nystrom approximation\n",
    "  def nystrom[T](x: Array[T], y: Array[Double], t: Array[T], kernel: MercerKernel[T], lambda: Double): KernelMachine[T]\n",
    "}\n",
    "```\n",
    "where the parameter `x` is the training data, `y` is the response variable, `kernel` is the Mercer kernel function, and `lambda` is the shrinkage/regularization parameter. The last two methods fit approximate models. For approximation models, the parameter t is the inducing input, which are pre-selected or inducing samples acting as active set of regressors. In simple case, these can be chosen randomly from the training set or as the centers of k-means clustering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cv.regression(10, x, y) { (x, y) =>\n",
    "  val t = smile.clustering.kmeans(x, 20).centroids\n",
    "  gpr.approx(x, y, t, new GaussianKernel(0.06), 0.01)\n",
    "}"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Scala (2.13)",
   "language": "scala",
   "name": "scala213"
  },
  "language_info": {
   "codemirror_mode": "text/x-scala",
   "file_extension": ".sc",
   "mimetype": "text/x-scala",
   "name": "scala",
   "nbconvert_exporter": "script",
   "version": "2.13.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
